function Mat = vdm23(d,b1,b2,b3)
% comput the Bform of degreed d at points (b1,b2,b3);
m = (d+1)*(d+2)/2;
plot_m = length(b1);
[I,J,K] = indices(d);
IM = diag(I)*ones(m,plot_m);
JM = diag(J)*ones(m,plot_m);
KM = diag(K)*ones(m,plot_m);

plot_IM = diag(b1)*ones(plot_m,m);
plot_JM = diag(b2)*ones(plot_m,m);
plot_KM = diag(b3)*ones(plot_m,m);
Mat = (plot_IM).^(IM').*(plot_JM).^(JM').*(plot_KM).^(KM');
IF = gamma(I+1);
JF = gamma(J+1);
KF = gamma(K+1);
A = factorial(d)*ones(plot_m,m)*diag(1./(IF.*JF.*KF));
Mat = A.*Mat;